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ABSTRACT 

A  new  method,  based  on  empirical  La  Place  transform, 
was  developed  to  find  asymptotically  efficient  estimates  of 
negative  binomial  distribution  parameters.   These  estimates 
were  found  fairly  close  to  those  found  by  the  method  of 
maximum  likelihood.   Efficiencies  over  95  percent  were 
obtained.   The  method  was  tested  with  a  set  of  data, 
generated  by  computer,  and  found  to  be  satisfactory  except 
in  a  few  cases.   Maximum  likelihood  also  fails  to  be 
satisfactory  in  these  cases. 
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I.   INTRODUCTION 

The  need  for  readily  computed  parameter  estimates  is 
great.   Maximum  likelihood  estimators  are  known  to  be  asymp- 
totically efficient,  but  in  many  settings  they  are  hard  to 
find.   The  most  popular  alternative  is  the  method  of  moments 
which  usually  yields  readily  computed  estimates,  but  some- 
times these  estimates  are  not  very  efficient.   See  ref .  [2] . 
This  study  looks  at  the  efficiency  of  a  method  which  uses 
the  probability  generating  function  (empirical  La  Place  trans- 
form evaluated  at  a  desirable  value  of  its  argument,  u  ). 

The  method  presented  herein  requires  computing  power 
greater  than  that  of  the  method  of  moments,  but  less  than 
that  of  maximum  likelihood,  which  calls  for  a  psi  function 
capability  .   It  can  be  used  with  a  hand  held  calculator  that 
has  the  lower  order  transcendental  functions,  i.e.,  logarithm, 
roots  and  power. 

The  basic  idea  is  to  select  a  system  of  estimating  equa- 
tions which  equates  various  statistics  to  their  expected 
values.   In  chapter  II,  estimating  equations  were  set  and 
efficiency  computations  were  done  using  theoretical  work 
presented  in  [3] .   Chapter  III  contains  the  procedure  to  be 
followed, applications  of  the  method  to  some  data  sets,  and 
comparisons  with  the  other  methods. 

In  chapter  IV,  efficiency  and  optimum  u   and  t  (=-log  u  ) 
tables  for  various  p-r  combinations  are  given  and  efficiency 
contours  in  three  different  planes  are  graphed. 


II.   METHODOLOGY 

A  negative  binomial  random  variable  X  has  the  probability 
function 


f(X-x;r,p)   =   ^.(r(+rf  qV  (2.D 


for  x  =  0,1,2,  ...   0<r,  0  <  p  <  1,  p  +  q  =  1.   The  partial 
derivatives  of  its  logarithm  are,  using  h  =  log  f 


9h     r  _  x 
3p      p    q 


||  =   T(r  +  x)  -  ?(r)  +  logp 


Using  the  basic  recursive  formula  for  the  psi  function  [1] 

x 
y(r  +  x)  -  ?(r)   =   I   r}1  +  i 

j=l 

one  may  then  express  the  system  of  maximum  likelihood 
equations  as 


x  -  sa   =    o 

p 


Ave      7  r— — r  +  log 

.  ,      *•  r  -  1  +  i     ^ 
i=l, . . . ,n        J 


(2.3) 
p   =   0 


where  the  values  x, ,  . . . ,  x   are  the  data  that  result  when 

1        n 

the  population  (2.1)  is  sampled  n  times.   The  system  (2.3) 
is  nonlinear  in  r,p  and  difficult  to  solve.   Iterative 
methods  must  be  used  and  the  second  member  of  (2.3)  (or 
some  version  thereof)  must  be  recomputed  in  each  cycle. 
This  is  the  main  computational  difficulty  in  using  maximum 
likelihood  in  this  setting. 

The  information  matrix  A,  which  is  defined  as: 


A   ■   -E 


afh 

ap2 


8r3p 


32h 

9p9r 


3r2 


is  needed  for  efficiency  calculations  and  can  be  developed 
using  the  methods  presented  in  [2] .   The  result  is 


1 

P 


where 


A   = 


qp 


l 
p 


l22 


A22   =   V  (r)  -  E(V  (r  +  x)  ) 


The  determinant  of  A  is  found  (using  [2] )  as 


1   V   q     .  r!nl  ,  (2  4) 

2  L    (n+i)    (r+n)  !  u,4; 


n=l 


The  two  estimating  equations  g,  and  g_,  which  are 


exploited  herein,  are  given  as 


g1(r,p;x)   =  x  "  p   " 


(2.5) 


g2(r,p;x)   =   tX  -  (j^)  *      =      0 


where 


n 
x. 


tX  =  i  I  t  * 


n 
i=l 


is  the  empirical  generating  function.   (Note:   the  substitution 
t  =  exp(-u)  converts  this  to  La  Place  transform).   Reduction 
to  a  single  equation  is  effected  by  solving  g,  =  0  for  p 
in  terms  of  x  and  r, 


P  =  zr= —  (2.6) 

x  +  r 


and  substituting  this  into  g~.   A  function  f(r)  is  obtained 
in  terms  of  known  variables,  namely  t  ,  x,  and  t, 
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f  (r)   =   ( -£ )r  -  tX   =   0       (2.7) 

r  +  x(l  -  t) 


We  can  solve  equation  (2.7)  for  r  by  Newton's  method,  the 
iterative  relationship  being 


f  (r±) 


ri+l   =   ri  "  f ' (r. ) 

l 


where 


f'(r)   =   ( £ )r  •  [(   x(^  t} )  +  log( £ )] 

r  +  x(l-t)        r  +  x(l-t)        r  +  x(l-t) 


assuming  an  initial  value  of  r  and  a  suitable  t  can  be 
found  (see  table  4.2).   After  solving  for  r,  p  can  be  found 
using  (2.6) . 

Efficiency  for  this  r-p  pair  is  given  by  (see  [2,3]) 


Eff     =     iii L  (2.8) 

|A| 


where 


-1  IaI2 

\lf±\       =      -lAL-  (2.9) 


|c 


|A|  is  given  in  (2.4),  the  matices  A  and  C  in  (2.9)  are 
defined  as 
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3g3 

3P 


3g: 


A   =   E 


3g. 


9g. 

IF 


C   =   nE 


gxg2 


gxg2 


and  g. ,g„  are  from  (2.5).   After  some  calculation,  A  and 
C  turn  out  to  be 


A      = 


-(-^-j^1    .    rd-t> 


C      = 


■1-qt 


=1 

P 


rqpr"1(l-t) 
d-qt)r+1 


d-qt) 


2 

P 


:rE-,r.l0g(i=at, 

1-qt  P 


(2.10) 


rqpr_1(l-t] 


(1-qt) 


r+1 


(-^T)r    -     ( 
1-qt2 


£_}2r 


1-qt 


(2.11) 


The  efficiency  of  the  estimation  scheme  is  found,  substituting 
|a| ,  |C|  and  lAl  into  (2.9)  and  (2.8),  as 
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BfiE-32 


ni^t,  lcg(i^  .  .aiizt),^ 


q     ,n   ,.2(r+l)         ~  „   «    n     .  , 

"    r  (1-qt)        r.,   ..2  r.,  ,.2,  ,  r    q     r!n!  -, 

[    ^2)r p  (1_gt)  _rqp  (1_t)  ][J17^TT--(iTK)T] 

(2.12) 


Detailed  computations  are  given  in  Appendix  A.   The  initial 
value  of  r  is  found  using  the  method  of  moments  (p,r), 
and  the  value  of  t  used  maximizes  the  efficiency  (2.12) 
when  evaluated  at  p,  r. 
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III.   USAGE  AND  APPLICATION 

In  this  chapter  the  procedure  to  be  followed,  in  order 
to  use  the  method  developed,  is  given  in  greater  detail  and 
is  applied  to  three  sets  of  data  taken  from  [4]  and  sixty 
sets  of  data  generated  by  computer  simulation  using  variqus 
combinations  of  the  parameters. 

A.   PROCEDURE 

The  steps  of  procedure  are  given  below. 

1.  Using  the  method  of  moments,  find  p  and  r  as 
starting  values.  Starting  values  are  given  by  (3.1) 
(see  [2]) 


p   =  x/s 


(3.1) 


r  = 


-2 
x 


s2-x 


—  2 

where  x  is  the  sample  mean  and  s   is  the  sample  variance. 

2.   Using  p,  r  found  in  step  1,  find  the  value  of  t, 
which  maximizes  efficiency  function  (2.12),  and  call  it 
t  .  and  also  find  the  efficiency,  which  is  maximum  with  this 
pair  p,  r  and  t  ,  and  call  it  EFF1 .   Tables  (4.1)  and  (4.2) 
have  been  prepared  for  this.   Single  variable  search  and 
golden  section  search  were  used  by  the  author  to  find  tQ 
and  EFF1. 


14 


3.  Using  r  as  a  starting  value  and  t  ,  found  in  step  2, 
find  r*  as  a  solution  of  (2.7),  and  p*  using  (2.6).   Newton's 
method  was  used  by  the  author  to  find  r*. 

4.  Step  2  can  be  repeated  with  p*  and  r*  replacing  p 
and  r,  and  using  a  new  t   (update  t   from  Table  (4.2)  if 
the  current  estimated  efficiency  (Table  4.1)  is  not  suffi- 
ciently high) .   This  step  involves  the  recomputations  of 

t   and  is  seldom  needed.   The  new  efficiency  using  p*,  r* 
and  updated  t   is  called  EFF2. 

5.  Stop  if  all  the  following  conditions  are  met. 


(1  -  EFF1)  <_  e1  ,    or    |EFF2  -  EFF1  |  <  e. 


and 


i*  - 


P   '  e-a' 


and 


r*  ~  rl  i  e4 


Then  p*  and  r*  are  the  estimators  of  true  parameters  p  and 
r,  and  the  estimated  efficiency  is  EFF2. 

6.   If  one  or  more  of  the  conditions  in  Step  5  are  not 
met,  let 


P  =   P* 
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EFF1   =   EFF2 


and  go  to  step  3  and  repeat  the  steps  following  until  the 
conditions  in  step  5  are  met.   e,  ,  e~,    e^  and  e.  are  the 
stopping  criteria  which  can  be  chosen  by  the  user . 


B.   EXAMPLES 

Examples  one  through  three  are  the  applications  of  the 
method  to  the  Cricket  score  data  of  Reep,  Polard  and 
Benjamin  [4] .   Example  four  is  based  on  data  generated  by 
computer  and  contains  sixty  cases. 

Example  1:   Applying  method  developed  to  the  Cowdrey 
data  [4]  the  following  information  is  obtained. 

Sample  mean  =   1.692 
Sample  variance  =   4.343 

p*   =   0.325 

r*   =   0.816 
Optimum  t    =   0.4  59 
Efficiency  =   0.996 
Number  of  Iterations   =   3 


/\  /\ 


With  the  same  data  maximum  likelihood  estimators  p,  r, 
and  method  of  moment  estimates  p,  r  are  given  in  [2]  as 
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p   =   0.329  r   -   0.831 

p   =   0.390  r   =   1.080 

Example  2:   Example  1  was  repeated  for  Barrington  data 
[4]  and  following  information  is  obtained. 

Sample  mean   =   2.095 

Sample  variance   =   4.939 

p*   =   0.346 

r*   =   1.111 

Optimum  t    =   0.538 
o 

Efficiency   =   0.995 

Number  of  Iterations  =   3 

p   =   0.345  r   =   1.014 

p   =   0.424  r   =   1.543 

Example  3 :   Like  the  previous  two  examples,  the  following 
information  is  obtained  by  applying  the  method  to  the 
Graveney  data  [4] : 

Sample  mean   =   1.57  0 
Sample  variance   =   4.474 

p*   =   0.315 

r*   =   0.722 
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Optimum  t    =   0.4  30 
o 

Efficiency   =   0.996 
Number  of  iterations   =   2 

p      =      0.317  r      =       0.729 

p      =      0.351  r      =      0.849 

When  the  first  three  examples  are  examined  carefully 
it  is  noticed  that  the  method  developed  in  this  study  is 
almost  as  efficient  as  maximum  likelihood  method.   Estimators 
are  pretty  close  to  maximum  likelihood  estimators  and 
much  better  than  those  found  by  the  method  of  moments. 

Example  4:   This  example  is  based  on  the  data  generated 
by  computer.   Sixty  cases,  each  of  which  has  a  different 
sample  size  and  parameters  p  and  r.   Four  sample  sizes, 
which  are  15,  30,  50  and  100;  three  r  values,  0.5,  2.5  and 
5;  and  five  p  values,  0.05,  0.1,  0.3,  0.5,  0.8  are  used. 
In  each  case  the  data  is  generated  for  given  sample  size 
(n) ,  p  and  r.   Results,  which  contain  case  number  (case  no.), 
sample  size  (n) ,  p,  r,  method  of  moments  estimates  (p,r), 
optimum  u  at  p-r  (Initial  u  ) ,  efficiency  at  the  end  of 
the  first  iteration  (Initial  eff.),  empirical  La  Place 
transform  estimates  (p*,r*),  optimum  u   at  p*-r*  (Final  u  ), 
efficiency  at  p*-r*-final  u   (final  eff)  and  number  of 
iterations  (No.  of  itr.),  are  tabulated  in  the  following 
pages.   Data  generated  for  the  sixty  cases  are  given  in 
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Appendix  B.   The  computer  program,  written  in  FORTRAN  IV, 
is  also  given  in  Appendix  C. 

When  all  the  cases  are  examined,  it's  seen  that  esti- 
mates are  not  given  for  cases  2,  3  and  18.   Case  3  provides 

nothing  to  work  with  as  all  samples  are  zero,  yielding 

—       2  2    — 

X  =  0,  s   =0.   The  other  two  cases  have  s   <  x  which 

signals  trouble  because  the  method  of  moments  estimates 

r  <  0,  p  >  1;  values  which  cannot  be  used  to  initialize 

our  iteration  scheme  (3.1)  or  the  maximum  likelihood  scheme 

(2.3)  . 

For  each  of  these  cases  the  function  f  (r)  of  (2.7) 

is  a  non-negative  decreasing  function  of  r  which  tends  to 

+0°  as  r  +   0  and  is  asymptotically  zero  as  r  ■*  ot.   Thus  no 

root  exists.   The  maximum  likelihood  approach  yields  a 

comparable  situation. 
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IV.   RESULT  AND  CONCLUSION 

In  this  chapter  efficiencies  for  various  p-r  pairs  were 
given  in  Table  (4.1).   The  method  developed  was  used  to 
calculate  table  entries.   Three  entries  were  given  for  each 
p-r  pair.   Entry  in  the  middle  is  efficiency  value  for  that 

pair  at  optimum  t   or  u  .   Entries  inside  parenthesis  are 

c  c  o      o  e 

also  efficiencies  and  found  as  follows: 

Let  u  =  -  log  t 
o       r   o 

If  we  miss  u  by  ±50  percent,  what  would  the  efficiency  be? 

Efficiency  inside  the  parenthesis  above  the  optimum  efficiency 

answers  the  question  for  u,  =  0.5  u  ,  and  other  one  below 
^  1        o 

the  optimum  answers  for  u-  =  1.5  u  . 
c  2        o 

The  change  u  =  -log  t  (or  rather  its  inverse  t  =  e   ) 
changes  the  probability  generating  t   function  into  the 
empirical  La  Place  transform.   It  is  seen  that  a  50  percent 
error  in  its  argument  u  results  in  only  a  minor  degradation 
of  efficiency,  whereas  a  50  percent  error  in  t  would  have 
disasterous  effects.   In  all  of  the  examples,  the  initial 
u  was  well  within  this  50  percent  range,  the  worse  case 
being  38  percent. 

After  examining  Table  (4.1)  one  can  conclude  that  effi- 
ciency is  monotonically  increasing  with  increasing  p  and 
increasing  r  for  r  greater  than  one  and  minimum  efficiency 

at  optimum  t  or  u   is  greater  than  0.95. 
c  o     o 
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Following  Table  (4.1)  optimum  t   and  u   (in  parenthesis) 
were  tabulated  in  Table  (4.2).   One  may  observe  that  t 
decreases  with  p  and  increases  with  r,  whereas  u   behaves 
oppositely. 

After  Table  (4.2)  efficiency  contours  were  graphed  in 
p-r  plane  (Figure  4.1),  in  mean-r  plane  (Figure  4.2)  and 
in  mean-exp(r)  plane  (Figure  4.3).   Compared  with  the 
Figure  in  [5]  Figure  (4.3)  always  gives  better  efficiency 
values  than  all  the  other  alternatives. 
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Figure  4.1   Efficiency  contours  in  p-r  plane 
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APPENDIX  A 

1 

DETAILED  COMPUTATIONS 

The  two  estimating  equations  g,  and  g2  given  in  chapter 
II  are 


9i  = 


x  -  SL  =   o 


*2  "  *   "  (I%>*   -   0 


Matrix  A's  entries  were  given  as  follows: 


'11 


8gl 


E(-%> 


12 


3gl 
3r 


E(-  S) 
P 


q 
p 


k21 


3g? 


"  =  <-<!%> 


p_.r-l  .  r(l-t) 


(1-qt) 


,    P  >r-l  r(l-t) 


l22 


E(Ir1) 


.tlnciat,^)'] 


^^)^)r 
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Matrix  C's  entries  were  given  as  follows 


'11 


n   E(gh      =      n   E[(x   -   ^)2] 

*^"  XT 


n 


=      n  Var(x)       =      n  Var(-      Y   x.) 

n      <*      1 

i=l 


n 


=     n--^-  Var(    \   x.)      =     Var(x) 
n 


i=l 


2 

P 


C12      =      n  E[gi.g2] 


nE[(x-f)(tx-    (rL,r,] 


n  Cov(X,tx) 


1      n  — 

n  Cov(-     I   xi7tx) 

i=l 


Cov(   xi,tx) 


n 


n 


n 


I   Cov(x. ,tX) 
i=l 


x  . 


I   Cov(x.,i      I    t   j) 
i=l  j=l 


n        n 


x  . 


i      I         I   Cov(Xi,t   h 
i=l   j=l 
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X  . 

Since  Cov(x.,t  -1 )  =  0  for  i  ^  j  because  of  independence, 
it  follows  that 


n 


x 


C    =   C 
c12      u21 


-    y  cov(x. ,t  x) 

n  L  1 

i=l 


■   Cov(x,tX)   =   E(x-tX)  -  E(x)E(tX) 


X 


One  can  calculate  E(x-t  )  using  the  following  relation. 


G(t)   =   E(tX)      G'  (t)   =   E(xtX-1) 


Then 


x 


E(xt  )   =   tG"  (t) 


where 


G(t)   ■   (i^t>r 


taking  the  derivative  and  multiplying  by  t,  E(xtx)  is  found 
as 


E(xtx)  =  i£3L_ 

d-qt) 


r+1 
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Then  it  yields 


r  =     r  =   rgpr   (1-t) 

12   "    21    "   (l-qt)r+1 


C22   =   n  E(g2)   .   E[(tX  -  (j^Vl 


n 
— —  -,      x. 

=  n  Var(t  )   =   n  Var(-  £  t  1) 

i=l 


n  * 

i=l 


Var(tx)   =  E(t2x)  -  [E(tX)]2 


=   G(t2)  -  [G(t)]2 


=   (   P   )r  -  (  P  )2r 
(l-qt2)      W 


As  the  reader  will  recall,  |A|  was  given  by  (2.4)  and  can  be 
written  as  follows. 


a|    =    ir  j  aL. r!nl     +    y     sC--£iaJL] 

1       2L  *  n+1   (n+r) !     *    n+1  (n+r) ! J 
p  n=l  n=ra+l 


Call  the  summation  from  n  =  m+1  to  °o  the  remainder.   The 
remainder  can  be  adjusted  such  that  it  will  be  smaller  than 
any  specified  small  value  (e) . 


45 


n     .  .  n 

y   3_  r!n!  <   y   q 

L        n+1  (n+r) !  -     L        n+1 
n=m+l  n=m+l 


r '  n ' 
Because  ,  '+r\t    ±   1  (reciprocal  of  binomial  coefficient)  , 


00  00 

Y        q q  y         n  _      q 

^         n+1  -     ~m+2      i      q  _       (m+2)p      -      e 
n=m+l                                  n=0 


A  value  of  m  can  be  found  such  that 

m+1 
(m+2)p  -  e 

holds.   Then  this  m  guarantees  the  remainder  will  be  smaller 
than  e,  and  the  efficiency  equation  of  (2.12)  becomes 


r+2                     [list  i„ (±=St,    -  SLiiztI]2 
Ef  f  =  £2 ,P  P     £ 


f(l-gt)2(r+1)        rM      ^>2  rn.,2w?qn        r!ni 

t  2   r p    (1_qt)    "rqp    (1_t)    ](    I   H+1    (r+n) ' 

(l+qtZ)  n=ln         ^         ; 
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APPENDIX  B 
SAMPLES  GENERATED  BY  COMPUTER 


Samples  generated  by  computer  for  sixty  cases  used  in 
chapter  three,  example  four,  are  given  in  the  following 
pages  case  by  case. 

See  subroutine  GNRNB  of  computer  program- 3  of  Appendix 
C  for  the  simulation  procedure,  which  was  taken  from  [6] . 
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VAL.  FRQ.  VAL.  FRQ.  VAL.         FRQ.  VAL.  FRQ. 


CASE  1 

CASE  7 

CASE 

10 

CASE 

13 

0       2 

3 

2 

(Cont'd) 

(Cont'd) 

1       7 

5 

1 

63 

1 

68 

1 

2       1 

8 

1 

64 

1 

70 

1 

3       3 

9 

1 

71 

1 

88 

1 

5       2 

10 

2 

73 

1 

120 

1 

11 

1 

74 

1 

131 

1 

CASE  2 

12 

13 

1 
1 

76 

1 

139 
153 

1 
1 

0       8 

14 
15 
19 
21 
22 

1 
1 
1 
1 
1 

2 

7 
8 

CASE 

11 

176 
12 

CASE 

1 

1  5 

2  2 

1 

1   • 
1 

14 

CASE  3 

1 

CASE  8 

9 
11 
17 

1 
1 
1 

15 
24 
30 

1 
1 

0      15 

1 

CASE  4 

1 
2 
3 

1 
2 
3 

18 
19 
22 

1 
1 
2 

40 
42 
50 

1 
1 

1 

0       1 

5 

1 

25 

1 

52 

1 

1       1 

6 

1 

27 

1 

59 

1 

2       1 

7 

1 

30 

1 

62 

1 

3       4 

8 

3 

37 

1 

70 

1 

4       1 

9 

1 

41 

1 

71 

1 

5       3 

11 

1 

78 

1 

7       2 

14 

1 

CASE 

12 

85 

1 

9       1 
12       1 

96 

1 
15 

CASE  9 

0 

4 

CASE 

CASE  5 

0 

1 

5 
2 

1 
4 
5 

2 

1 
1 

0 

3 

0       2 

2 

2 

6 

3 

1 

2 

1       1 

4 

2 

7 

1 

2 

1 

2       4 

5 

2 

10 

1 

3 

1 

3       1 

6 

2 

11 

1 

5 

1 

4       2 

25 

1 

9 

2 

5       1 

CASE  10 

15 

1 

7       4 

CASE 

13 

18 
22 
30 
34 

1 
1 
1 
1 

11 
14 
21 

1 
1 
1 

33 
49 

CASE  6 

1 
1 

0  8 

1  5 

30 
33 
39 

1 
2 

1 

52 
59 
60 

1 
1 
1 

CASE 

16 

2       1 

5       1 

40 

1 

63 

1 

0 

8 

57 

1 

64 

1 

1 
2 
3 
4 
5 

7 
6 
6 
1 
2 
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VAL.    I 

?RQ. 
17 

VAL.    FRQ. 
CASE  22 

VAL 

FRQ. 

VAL 

FRQ. 

CASE 

14 

CASE  25 

0 

CASE  27 

0 

14 

1 

1 

1 

8 

1 

8 

3 

2 

23 

1 

1 

4 

2 

5 

6 

1 

28 

3 

2 

1 

3 

3 

7 

3 

29 

1 

3 

3 

8 

4 

30 

1 

4 

3 

CASE 

18 

9 

2 

31 

1 

5 

2 

11 
12 

5 

1 

32 
33 

1 
4 

7 
8 

1 
1 

0 

1 

27 
3 

13 

1 

34 

1 

9 

2 

14 

2 

36 

1 

11 

1 

CASE 

19 

15 
16 
17 

1 
1 
1 

37 
47 
51 

2 

1 
1 

15 
18 
26 

1 
1 

1 

1 

2 

19 

2 

53 

1 

29 

1 

2 

6 

20 

2 

58 

1 

3 

7 

26 

1 

60 

1 

CASE  2  8 

4 

5 

6 

7 
10 
11 
14 

CASE 

3 
4 
3 

1 
1 
2 

1 

20 

67 
68 
71 
73 
77 
81 
88 
102 

1 
1 
1 
1 
1 
1 
1 
1 

CASE  26 

CASE 

23 

42 
46 
52 
55 
57 
58 
59 
60 
65 
72 

1 

0 
1 
2 
3 
4 
5 
6 
7 

3 
2 
2 
3 

4 
5 
4 

1 

1 

1 
1 
1 
1 
1 
1 

5 

1 

0 

1 

1 

8 

8 

2 

0 

1 

73 

1 

2 

2 

10 

2 

5 

1 

79 

1 

3 

2 

15 

1 

6 

1 

81 

1 

4 

7 

20 

1 

8 

1 

91 

1 

5 

4 

9 

1 

92 

1 

7 

2 

CASE 

24 

10 

1 

98 

1 

12 
15 
16 
17 
20 
21 

1 
3 
2 

1 
2 

1 

100 
106 
112 
118 
120 
121 

1 
1 
1 
1 
1 
1 

CASE 

21 

0 

12 

0 
1 
2 
3 
4 

19 
5 
3 
2 

1 

1 
2 
3 
4 
5 
7 

10 
2 

1 
1 
2 
2 

23 

1 

124 

1 

24 

1 

129 

1 

27 

2 

132 

1 

29 

1 

141 

1 

31 

1 

148 

1 

32 

1 

167 

1 

33 

1 

185 

1 

35 

1 

197 

1 

47 

1 

49 

1 

53 

1 

61 

1 
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VAL.  FRQ.  VAL.  FRQ.  VAL.  FRQ.  VAL.  FRQ. 


CASE 

29 

0 

CASE 

31 

CASE 

36 

CASE 

39 

7 

1 

21 

0 

37 

0 

27 

12 

1 

1 

13 

1 

6 

1 

5 

15 

1 

2 

11 

2 

6 

2 

8 

18 

1 

3 

4 

3 

1 

3 

4 

21 

3 

7 

1 

4 

3 

22 

1 

CASE 

37 

5 

1 

26 

1 

CASE 

32 

6 

1 

30 

1 

2 

1 

10 

1 

33 

40 
49 
52 
54 
55 
57 
58 
65 

1 

1 
1 
1 
4 
2 
1 
1 
1 
1 

0 
1 
2 
3 
4 

24 

17 

6 

1 
1 

3 
4 
5 
6 
7 

3 
3 
2 
2 
2 

CASE 

40 

8 
16 

1 
1 

5 

CASE 

1 
33 

8 

9 

10 

11 

12 

5 
3 
3 
2 

6 

17 
22 
26 
27 
28 

1 
1 

3 

4 

66 

1 

68 
71 

1 
2 

0 

1 

44 
4 

13 
14 

1 
2 

30 

31 

1 
1 

74 
86 

1 

1 

2 

2 

15 
17 

4 
4 

32 
33 

3 
3 

95 

1 

CASE 

34 

20 
22 

1 
2 

35 
36 

2 

2 

CASE 

30 

1 
2 

8 
6 

23 
25 

1 
1 

37 
40 

1 

2 

0 

6 

3 

3 

29 

1 

41 

1 

1 

5 

4 

7 

33 

1 

44 

2 

2 

2 

5 

10 

45 

1 

3 

2 

6 

4 

CASE 

38 

47 

1 

4 
5 

1 
2 

7 
8 

2 

4 

48 
49 

1 
1 

0 

1 

2 

3 

4 

5 

6 

7 

8 

10 

11 

14 

1 
7 
4 
3 
5 
6 
4 
6 
2 
3 
2 
1 

6 

1 

9 

4 

54 

2 

8 

1 

11 

1 

55 

1 

9 

1 

17 

1 

58 

1 

10 

1 

60 

1 

12 

2 

CASE 

35 

61 

1 

15 

1 

62 

1 

21 
25 
38 
47 
52 

1 
1 
1 
1 

1 

0 
1 
2 
3 
4 

9 
9 
8 
8 

5 

68 
70 
71 
74 
75 
87 

1 
1 
2 
1 
1 
1 

5 

6 

15 

2 

91 

1 

6 

1 

16 

2 

93 

1 

7 

2 

17 

1 

8 

1 

25 

1 

9 

1 

50 


VAL. 


FRQ. 


VAL. 


FRQ. 


VAL. 


FRQ. 


VAL, 


FRQ. 


CASE 

41 

CASE 

43 

CASE 

44 

CASE 

45 

(Cont'd) 

3 

1 

17 

1 

7 

1 

4 

2 

21 

1 

8 

2 

8 

3 

8 

3 

33 

1 

12 

1 

9 

1 

9 

3 

40 

1 

13 

1 

10 

4 

10 

3 

43 

1 

15 

1 

11 

3 

11 

1 

.  45 

1 

21 

3 

12 

1 

12 

1 

50 

1 

22 

1 

13 

1 

13 

1 

51 

2 

24 

2 

14 

1 

14 

4 

52 

1 

25 

2 

23 

2 

15 

4 

55 

1 

27 

1 

26 

1 

16 

1 

57 

1 

29 

3 

31 

1 

17 

2 

61 

1 

30 

2 

34 

1 

18 

3 

63 

2 

32 

1 

38 

1 

19 

1 

64 

1 

35 

2 

59 

1 

20 

2 

68 

1 

38 

1 

64 

1 

21 

1 

69 

1 

40 

1 

22 

1 

71 

1 

41 

1 

CASE 

46 

23 

1 

76 

1 

42 

1 

24 

1 

78 

1 

43 

1 

0 
1 
2 
3 

35 

26 

1 

80 

1 

44 

1 

32 

27 

3 

82 

1 

45 

1 

20 

29 

1 

83 

1 

46 

1 

7 

33 

1 

87 

1 

47 

1 

4 
5 

5 

35 

1 

90 

1 

48 

1 

•J 

1 

41 

1 

93 

1 

49 

1 

47 

1 

94 

1 

51 

2 

CASE 

47 

50 

2 

1 

95 
98 

1 
1 

53 

55 

1 
1 

74 

76 

1 

100 

1 

56 

1 

0 

59 

82 

1 

101 

1 

70 

1 

1 

27 

102 

1 

73 

1 

2 

11 

CASE 

42 

104 

1 

77 

1 

3 

_2 

109 
114 

1 
2 

82 
85 

1 
1 

4 

1 

0 

1 

16 
4 
9 
2 

T! 

121 

1 

86 

1 

CASE 

48 

125 

1 

105 

1 

2 
3 

127 

1 

111 

1 

0 

92 

4 
5 

1 
3 

129 

1 

130 

1 

1 
2 

7 

135 

1 

140 

2 

1 

6 
7 

2 
1 
3 
2 

140 

1 

149 

1 

CASE 

45 

158 

1 

8 
10 

163 

1 

0 

12 

12 

3 

167 

1 

1 

4 

16 

2 

171 

1 

2 

3 

20 

1 

184 

1 

3 

5 

35 

1 

185 

1 

4 
6 
7 

2 

1 
1 
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VAL.  FRQ.  VAL.         FRQ.  VAL.  FRQ.  VAL.  FRQ. 


CASE 

49 

CASE  5  2 

CASE 

54 

CASE  55 

(Cont'd) 

(Cont'd) 

0 

1 

0 

60 

1 

8 

10 

7 

1 

15 

56 

1 

2 

12 

11 

2 

2 

10 

59 

5 

3 

13 

12 

5 

3 

4 

61 

1 

4 

13 

13 

5 

4 

6 

62 

2 

5 

17 

14 

4 

5 

2 

64 

1 

6 

11 

15 

5 

6 

1 

65 

1 

7 

9 

16 

5 

8 

1 

68 

3 

8 

4 

17 

4 

14 

1 

72 

4 

9 

2 

18 

1 

74 

1 

10 

2 

19 

1 

CASE 

55 

76 

1 

11 
15 

5 
2 

20 
21 

9 

1 

77 
78 

1 

13 
15 
18 
19 
22 
23 
24 
25 

1 
3 
1 
1 
1 
1 
2 
4 

J- 

1 

20 

1 

22 

2 

81 

2 

CASE 

50 

24 
26 

2 

1 

88 
90 

1 
1 

27 

1 

91 

1 

0 

15 

29 

1 

97 

2 

1 

21 

101 

1 

2 

22 

CASE  53 

112 

1 

3 

17 

118 

1 

•J 

4 
5 

12 
5 

0 

5 

26 
28 

2 

1 

CASE  56 

6 

7 

8 

13 

2 
2 
3 

1 

1 
2 

9 
3 

29 
30 

1 

1 

3 

10 

31 

1 

1 

1 

4 

16 

32 

2 

4 

2 

5 

9 

33 

2 

5 

7 

CASE 

51 

6 

7 

34 

1 

6 

1 

7 
8 

7 
7 

35 
36 

3 

6 

7 
8 

2 

3 

0 

70 

9 

3 

37 

1 

9 

2 

1 

18 

10 

7 

38 

4 

10 

4 

2 

8 

11 

5 

39 

3 

11 

2 

3 

1 

12 

2 

40 

2 

12 

3 

4 

2 

13 

2 

41 

1 

13 

2 

6 

1 

14 

2 

42 

2 

14 

4 

15 

2 

43 

1 

15 

2 

CASE 

52 

17 

1 

44 

2 

16 

4 

19 
23 

1 
1 

45 
46 

1 
3 

17 
18 

3 
2 

0 
2 
3 
4 
5 
6 
7 
8 
9 

1 
3 
3 
4 
4 
3 
7 
12 
7 

24 

1 

47 

2 

19 

2 

49 

1 

20 

4 

50 

2 

21 

1 

51 

1 

22 

3 

52 

1 

23 

4 

53 

3 

24 

5 

54 

2 

25 

2 

55 

2 

26 

1 

52 


VAL.  FRQ.  VAL.         FRQ.  VAL.  FRQ.  VAL.  FRQ. 


CASE  56 

CASE 

58 

CASE  58 

CASE  59 

(Cont'd) 

(Cont'd) 

(Cont'd) 

21 

1 

27 

4 

33 

2 

116 

1 

30 

3 

28 

2 

35 

1 

118 

2 

32 

3 

29 

2 

37 

1 

119 

1 

34 

2 

30 

2 

43 

1 

121 

1 

35 

5 

31 

3 

44 

2 

122 

1 

36 

1 

33 

3 

45 

4 

124 

1 

37 

2 

34 

3 

46 

4 

125 

1 

38 

1 

35 

1 

50 

1 

127 

2 

39 

1 

36 

1 

53 

1 

131 

1 

40 

2 

37 

1 

54 

1 

135 

1 

41 

2 

38 

1 

55 

2 

137 

1 

42 

2 

39 

2 

57 

1 

140 

1 

44 

1 

40 

2 

60 

1 

141 

1 

45 

1 

46 

1 

64 

1 

144 

1 

47 

1 

48 

1 

65 

1 

147 

1 

49 

2 

50 

1 

68 

1 

152 

1 

51 

1 

53 

1 

69 

3 

153 

1 

52 

1 

54 

1 

71 

1 

161 

1 

53 

1 

72 

1 

72 

2 

180 

1 

54 

1 

76 

1 

73 

1 

185 

2 

56 

1 

74 

2 

190 

1 

57 

2 

CASE  57 

76 

2 

192 

1 

58 

1 

77 
78 

2 

1 

193 
197 

1 
1 

61 
63 

1 
1 

0 
1 
2 

32 
10 
10 
6 
9 
7 
2 
2 
4 
1 
2 
2 
1 
1 
3 
3 
1 
1 
1 
1 
1 

79 
80 

1 
1 

CASE  59 

64 
65 

1 
1 

81 

1 

66 

1 

3 

4 

5 

6 

7 

8 

9 

10 

12 

13 

14 

17 

18 

21 

23 

27 

29 

35 

82 

83 

85 

88 

90 

91 

93 

95 

96 

97 

99 

100 

101 

103 

104 

105 

108 

109 

1 
1 
2 
2 
2 
1 
1 
1 
2 
1 
1 
3 
3 
1 
1 
1 
2 
2 

7 
10 
11 
12 
13 
15 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 

2 

1 
2 
1 
1 
1 
1 
1 
1 
2 
2 
4 
3 
2 
2 
1 
2 
1 
3 

67 

69 

70 

72 

73 

78 

81 

82 

85 

88 

90 

91 

93 

98 

102 

103 

107 

110 

2 

1 
1 
1 
1 
2 
1 
3 
2 
1 
1 
1 
2 
1 
1 
1 
1 
1 

■ 

110 
114 

1 
1 

140 

1 

53 


VAL .  FRQ . 


CASE 

60 

0 

18 

1 

10 

2 

7 

3 

6 

4 

5 

5 

4 

6 

1 

7 

3 

8 

2 

9 

4 

10 

1 

11 

5 

12 

1 

13 

1 

14 

2 

15 

1 

16 

4 

17 

1 

18 

2 

19 

2 

20 

2 

21 

2 

22 

2 

25 

1 

26 

3 

28 

1 

30 

2 

31 

2 

32 

1 

38 

1 

50 

1 

53 

1 

90 

1 

54 


APPENDIX  C 
COMPUTER  PROGRAMS 

Computer  programs,  to  prepare  Table  (4.1)  and  Table  (4.2) 
(Computer  program  1) ,  to  find  the  estimates  for  the  examples 
one  through  three  of  chatper  three  (computer  program  2) , 
and  to  generate  the  data  for  the  sixty  cases  of  the  example 
four  of  chapter  three  and  to  find  the  estimates  for  them 
(computer  program  3) ,  are  given  in  the  following  pages. 
There  are  two  subroutines,  namely  GOLDEN  (Golden  section 
search) ,  SVS  (Single  variable  search) ,  and  one  function, 
namely  FCN  (to  evaluate  efficiency  function) ,  which  are 
common  to  the  three  programs  mentioned  above.   These  are 
given  at  the  beginning  under  the  title  COMMON  SUBROUTINES. 
After  that  the  three  programs  are  given. 

One  who  wants  to  use  any  one  of  the  programs,  must  include 
the  common  subroutines  after  the  main  program  of  the  computer 
program  he  chooses . 
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COMMON  SUBROUTINES 


SUBROUTINE  SVS ( FCN , TNOT ) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

COMMON  P,R,DETLAM 

Q=1.0DO-P 

UPPER=1 . ODO/DABS (Q) 

ALOW=-1.0DO/DSQRT(DABS (Q) ) 

Riou=o.oooooinn 

START=0.9999 

U= START 

S=0.0001D0 

A=FCN(U) 

U1=U+S 

B=FCN(U1) 

IF  (B.GT.A)  GO  TO  5 

IF  (B.EQ.A)  GO  TO  6 

M=3 

1  U2=START+S*2.0D0** (M-2) 
IF  (U2.GE. UPPER)  GO  TO  99 
A=B 

B=FCN(U2) 

IF(A.EQ.B)  GO  TO  7 

IF  (A.LT.B)  GO  TO  8 

M=M+1 

U=U1 

U1=U2 

GO  TO  1 

7  CALL  GOLDEN (FCN, Ul,U2,RIOU, 50, TNOT) 
GO  TO  999 

99     CALL  GOLDEN ( FCN, Ul, UPPER, RIOU, 50, TNOT) 
GO  TO  999 

8  CALL  GOLDEN ( FCN, U,U2, RIOU, 50, TNOT) 
GO  TO  999 

6      M=2 

CALL  GOLDEN (FCN , U , Ul , RIOU ,50, TNOT ) 

GO  TO  999 
5      M=3 

2  U2=START-S*2.0D0** (M-3) 
IF  (U2.LE.ALOW)  GO  TO  7  7 
B=A 

A=FCN(U2) 

IF  (A.EQ.B)   GO  TO  17 

IF  (A.GT.B)   GO  TO  18 

U1=U 

U=U2 

M=M+1 

GO  TO  2 

17  CALL  GOLDEN (FCN, U2,U, RIOU, 50, TNOT) 
GO  TO  999 

18  CALL  GOLDEN (FCN, U2,U1, RIOU, 50, TNOT) 
GO  TO  999 

77     CALL  GOLDEN (FCN, ALOW, U, RIOU, 50, TNOT) 
999    RETURN 
END 
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SUBROUTINE  GOLDEN ( FCN , ENDL , ENDR , RIOU , NMAX , XOLD ) 
IMPLICIT  REAL* 8  (A-H,0-Z) 
COMMON  P,R,DETLAM 

SIGMA=(3.0D0-DSQRT(5.ODO) )/2.0D0 
AL=ENDR-ENDL 
-—  XOLD=ENDL+SIGMA*AL 
FOLD=FCN(XOLD) 
DO  1  1=2, NMAX 
XNEW=ENDL+ENDR-XOLD 
FNEW=FCN(XNEW) 
IF  (FNEW.LE.FOLD)  GO  TO  5 
IF  (XNEW.LT.XOLD)  GO  TO  6 
ENDR=XNEW 
20     AL=ENDR-ENDL 
GO  TO  30 

6  ENDL=XNEW 
GO  TO  20 

5      IF  (XNEW.LT.XOLD)  GO  TO  7 

ENDL=XOLD 
8      XOLD=XNEW 

FOLD=FNEW 

GO  TO  20 

7  ENDR=XOLD 
GO  TO  8 

30     IF(AL.LE.RIOU)  GO  TO  40 
1      CONTINUE 
4  0     RETURN 
END 


FUNCTION  FCN(X) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

COMMON  P,R,DETLAM 

ZR=1.0D-73 

Q=1.0D0-P 

A=R*P** (R+2.DO)/Q 
U-  B=1.0D0-Q*X 

C=(DABS ( (B*DLOG(B/P)/P-Q* (1.0D0-X)/P) ) )**2.D0 

D=B** (2.D0*R+2.0D0)/(1.0D0-Q*DABS(X) **2.D0)**R 

E=P**R*B**2.D0+R*Q*P**R* (DABS (1.0D0-X) ) **2.D0 

G=-A*C/DETLAM 

F=D-E 

IF  (DABS(F) .LT.ZR)  GO  TO  5 

FCN=G/F 

GO  TO  7 
5      FCN=-1.0D0 
7      RETURN 

END 


57 


COMPUTER  PROGRAM  1 


IMPLICIT  REAL*8  (A-H/0-Z) 

EXTERNAL  FCN 

DIMENSION  PR(14) ,RR(26) ,EFF(26,14) ,TZERO(26,14) ,TTEST(52) ,EFF1(26, 
114),EFF2(26,14) 

COMMON  P,R,DETLAM 

DATA  EFF,TZERO/728*0.0D0/ 

DATA  EFF1,EFF2/728*0.0D0/ 

DATA  PR/.05D0, .06D0, .07D0, .08D0, .09D0, .1D0, .2D0, .3D0, .4D0, .5D0, .6D 
10,.7D0,.8D0,.9D0/ 

DATA  RR/.5D0,1.D0,1.5D0,2.D0,2.5D0,3.D0,3.5D0,4.D0,4.5D0,5.D0,6.D0 
1,7.00,8.00,9.00,10.00,11.00,12.00,13.00,14.00,15.00,16.00,17.00,18 
l.D0,19.D0,0.0D0,0.0D0/ 

DATA  TTEST/52*0.0D0/ 

WRITE    (6,444) 
444         FORMATCl') 

EPS=0.0001D0 

DO  1  J=l,14 

P=PR(J) 

Q=1.0D0-P 

K=l 

7  TEST=Q**(K+1)/CP*(DFLQAT(K)+2.0D0)) 
K=K+1 

GO  TO  7 

8  M=K 

WRITE  (6,21)  M 
21     F0RMAT(15X, 'UPPER  BOUND=' ,18) 

DO  2  1=1,24 

R=RR(I) 

DETIAM=O.ODO 

DO  10  N=1,M 

RATTONKL.ODO 

DO  12  JJ=1,N 
12   '         RATia^RATION*DFL0AT(JJ)/(DFLGAT(JJ)+R) 

DETIAM=DETLAM40^*N*RATI0N/DFLC!AT  (N+l) 
10     CONTINUE 

CALL  SVS(FCN,TNOT) 

TZERO(I,J)=TNOT 

TPL=DEXP ( 0 . 5D0*DLOG (TNOT) ) 

TMIN=DEXP (1 . 5D0*DLOG (TNOT) ) 

EFF2 (I, J)=-FCN(TPL) 

EFF1 (I ,  J) =-FCN (TMIN) 

EFF (I , J) =-FCN (TNOT) 
2      CONTINUE 
1      CONTINUE 

WRITE  (6,200)   (PR(I) ,1=1,14) 

200    FORMAT (///,50X,'P  VALUES' ,/,12X,14 (2X,F5. 3, IX) ,/,12X,14 (2X, ' 

1\1X),/) 

DO  13  1=1,24 

WRITE  (6,202)  (EFF1(I,J) ,J=1,14) 

WRITE  (6,201)  RR(I),(EFF(I,J),J=1,14) 

WRITE  (6,202)  (EFF2(I,J) ,J=1,14) 
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201  FORMAT(3X, 'R=' ,F5.2,2X,14 (2X,F5.3,1X) ) 
WRITE  (6,203)  (TZERO(I,J) ,J=1,14) 

202  FORMAT  (12X,  14  (IX,  '  (',F5.3,')  ')) 

203  FORMAT(5X, 'OP.TNOT' ,14 (2X,F5.3,1X) ,//) 
13  CONTINUE 

STOP 
END 
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COMPUTER   PROGRAM      2 

IMPLICIT  REAL*8(A-H,0-Z) 
EXTERNAL  FCN 

DIMENSION  FRQ(20) , OBS (20) 
COMMON  P,R,DETLAM 
WRITE    (6,444) 
444  FORMAT('l') 

EPS=0.0001D0 
WRITE    (6,91) 

91  FORMAT (5X, 'READ  IN  NO.   OF  OBS.  AND  VARIATES  IN  215  FORMAT1,/) 
READ    (5,81)   NOBS,.NVAR 

81  FORMAT (215) 
WRITE    (6,92) 

92  FORMAT (5X,  'NOW  GIVE  ME  VARIATE  VALUES  AND  CORRESPONDING  FREQUENCY 
UN  THE  FORMAT  2F10.5  FOR  EACH  CARD' ,/) 

SUM=0.0D0 

SUM1=0.0D0 

DO  4  I=1,NVAR 

READ    (5,82)    CBS(I)  ,FRQ(I) 

82  FORMAT (2F10. 5) 
SUM=SUM+OBS (I) *FRQ(I) 
SUM1=SUM1+ (OBS (I) **2 . 0D0) *FRQ  (I) 
XBAR=SUM/DFLOAT (NOBS) 

SSQR=  (SUMl-DFLOAT (NOBS) *XBAR**2 . 0D0) /DFLOAT (NOBS-1) 
WRITE    (6,209)    XBAR,   SSQR 
209         FORMAT (15X, 'USING  GIVEN  DATA  FOLLOWING  INFORMATION  IS  OBTAINED ' ,// 
1,20X, 'SAMPLE  MEAN  =' ,F7. 4,//, 2 OX, 'SAMPLE  VARIANCE     =',F7.4,// 

2) 
P=XBAR/SSQR 

R=XBAR**2 . 0D0/ (SSQR-XBAR) 
Q=1.0D0-P 
K=l 

7  TEST=Q** (K+l) / (DFLOAT (K+2) ) 
IF  (TEST.LE.EPS)    GO  TO  8 
K=K+1 

GO  TO  7 

8  M=K 
DETLAM=0.0D0 
DO  10  N=1,M 
RATION=1.0D0 
DO  12  JJ=1,N 

12     RATION=RATICN*DFLOAT(JJ)/ (DFLOAT (JJ)+R) 

DETLAM=DETLAM4Q**N*RATION/D  FLOAT  (N+l ) 
10     CONTINUE 

L=l 
95     CALL  SVS (FCN,TNOT) 

EFF1=-FCN (TNOT) 

RSTAR=R 

PSTAR=P 

SUM2=0.0D0 

DO  14  J=1,NVAR 
14     SUM2=SUM2+FRQ  (J) *TNOT**OBS (J) 

TXBAR=SUM2/DFLOAT (NOBS) 
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CALL  CQMPR(TNGT,TXBAR,XBAR,RSTAR) 

P=R/  (R+XBAR) 

Q=1.0D0-P 

K=l 
77     TEST=Q**  (K+l)  /  (£*DFLOAT  (K+2) ) 

IF  (TEST.LE.EPS)  GO  TO  88 

K=K+1 

GO  TO  77 
88     M=K 

DETLAM=0.0D0 

DO  51  N=1,M 

RATTON=1.0D0 

DO  52  JJ=1,N 
52  RATIQ^I^TICN*DFLQAT(JJ)/(DFLQAT(JJ)+R) 

51  DFHIAl^DFJIi^H<2**N*RATION/DFI£lAT(N+l) 

EFF2=-FCN (TNOT) 

IF (DABS (P-PSTAR) .GT. 0.005)    GO  TO  29 

IF (DABS (R-RSTAR) .GT. 0.005)    GO  TO  29 

IF (DABS (EFF1-EFF2) .GT. 0.04)    GO  TO  29 

WRITE    (6,94)    P,R,EEF2,TNOT,L 

FORMAT(20X, 'ESTIMATE  OF  P         =' ,F7. 5,//, 20X, 'ESTIMATE  OF  R         =',F 
17. 4,//, 2 OX, EFFICIENCY  =' ,F7.5,//,20X, 'OPTIMUM  TNOT  =',F 

27.5,//,20X, 'NO.OF  ITERATIONS  =',FI4,//) 

GO  TO  999 
29  L=L+1 

GO  TO  95 
999         STOP 

END 


SUBROUTINE  CQMPR(T,TX,XB,RST) 
IMPLICIT  REAL*8    (A-H,0-Z) 
COMMON  P,R,DETLAM 
R1=RST 
30  B+R1+XB* (1 . 0DO-T) 

TY=^  (Rl/B)  **R1 
F=TX-TY 

FP1=-TY 

FP2=  (XB*  (1 .  0D0-T)  /B)  +DL0G  (Rl/B) 

FP=FP1*FP2 

R2=R1-F/FP 
IF(DABS(R2-R1).LE.0.00001D0)   GO  TO  66 
R1=R2 
GO  TO  30 
66  R=R2 

RETURN 
END 
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COMPUTER   PROGRAM      3 

IMPLICIT  REAL*8    (A-H,0-Z) 
EXTERNAL  FCN 

DIMENSION  FRQ(200)/OBS(200),PVECT(10) ,RVECT(10) ,NVECT(10) ,DAT(200) 
1,ISD(200) 
COMMON  P,R,DETLAM 

DATA  RVECT/5 . 0D0 , 2 . 5D0 , 0 . 5D0 , 7*0 . 0D0/ 
DATA  PVECT/.8D0,.5D0,.3D0,.1D0,.05D0,5*0.0D0/ 
DATA  NVECT/15,30,50,100,6*0/ 
CALL  OVFLOW 
WRITE (5, 541) 
541- .      FORMAT (//,15X, 'READ  ISEED  USED  IN  GENERATING  DATA') 
READ (5, 978)      IS 

978  FORMAT  (110) 
WRITE (6, 979)      IS 

979  FORMAT (15X, 'ISEED  USED  TO  GENERATE  DATA=', 110) 
CALL  INT  (13,130,200) 

WRITE    (6,444) 
444  FORMAT('l,) 

EPS=0.0001D0 

NCASE=0 

DO  112  IN=1,4 

DO  113  IP=1,5 

DO  114  IR=1,3 

P=PVECT(IP) 

R=RVECT(IR) 

NOBS=NVECT (IN) 

NCASE=NCASE+1 

NCP=2*NCASE 

NCM=NCP-1 

ISEED1=ISD (NCP) 

ISEED2=ISD (NCM) 

CALL  GNRNB(NOBS,OBS,P,R,ISEED2,ISDEEDl) 

OBMX=-99999.D0 

DO  321  LK=l,NOBS 

IF    (OBS(LK)  .GT.OBMX)      OBMX=OBS (LK) 
321         CONTINUE 

OBLMOBMX+IO.DO 

NVAR=0 
390         OBMN=OBMX+5.D0 

DO  370  MI=l,NOBS 

IF  (OBS (MI) . LT . OBMN)   OBMN=0BS  (MI ) 
370    CONTINUE 

ICNT=0 

DO  380  MJ=l,NOBS 

IF  (OBS (MJ).NE. OBMN)   GO  TO  380 

ICNT=ICNT+1 

OBS(MJ)=OBLM 
380    CONTINUE 

NVAR=NVAR+1 

DAT(NVAR)=OBMN 

FRQ  (NVAR) =DFLQAT (ICNT) 

IF(OBMN.EQ.OBMX)   GO  TO  395 
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GO  TO  390 
395    WRITE  (6,889) 

889    FORMAT (//,10X, 89 ('*'),/,40X, 'DATA  GENERATED  IS  GIVEN  BELOW' ,/,39X, 
131C*'),/) 
WRITE  (6,301)   (DAT(L0),L0=1,NVAR) 

301  FORMAT (20X, 'VALUES', /,20X, 6 ('-'), /, (5X,20(1X,F5.0) ,/) ) 
WRITE  (6,302)   (FRQ(LO) ,L0=1,NVAR) 

302  fonrat920X, 'FREQUENCIES', /,20X,11( '-' ) ,/, (5X,20(1X,F5.0) ,/) ) 
WRITE(6,543) 

543  FORMAT (//,10X,' DATA  WRITTEN  ABOVE  IS  BASED  ON  FOLLOWING  PARAMETERS 
l',/,10X,  52  (■-•)■) 

WRITE  (6,544)  NCASE,P,R,NOBS,NVAR 

544  FORMAT (/,10X, 'CASE  NO.=' ,15. 5X, 'P=' ,F7.4,5X, 'R=' ,F7.2,5X, 'N=' ,Ir,5 
IX, 'NO. OF  VALUES=' ,13,/) 

SUM=0.0D0 
SUM1=0.0D0 
DO  445  I=1,NVAR 
SUM=SUM+FRQ  (I) *DAT (I) 
445    SUM1=SUM1+FRQ(I)*DAT(I)**2.0D0 
XBAR=SUM/DFLOAT (NOBS ) 

SSQR=  (SUMl-DFLOAT (NOBS) *XBAR**2 . 0D0) /DFLOAT (NOBS-1) 
IF  (XBAR.GE.SSQR)  GO  TO  71 
P=XBAR/SSQR 

R=XBAR**2 . 0D0/ (SSQR-XBAR) 
GO  TO  72 

71  R=5.0D0 
P=0.9D0 
GO  TO  114 

72  Q=1.0D0-P 
K=l 

7  TEST=Q**(K+1)/ (DFLOAT (K+2) ) 
IF (TEST.LE.EPS)  GO  TO  8 
K=K+1 

GO  TO  7 

8  M=K 
DETLAM+0.0D0 
DO  10  N=1,M 
RATION=1.0D0 
DO  12  JJ=1,N 

12  RATION=RATICN*DFLOAT(JJ)/ (DFLOAT (JJ)+R) 

DETILAM=DE[IAi:4+<2**N*RATiasI/rFLOAT  (N+l) 
10  CONTINUE 

L=l 
95  CALL  SVS (FCN/TNOT) 

EFF1=-FCN (TNQT) 

RSTAR=R 

PSTAR=P 

SUM2=0.0D0 

DO  145  J=1,NVAR 
145  SUM2=SUM2+FRQ(J) (TNOT**DAT ( J) 

TXBAR=SUM2/DFLQAT (NOBS) 

CALL  COMPR('I,NCT/TXBAR,XBAR,RSTAR) 
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P=R/  (R+XBAR) 
Q=1.0DO-P 
K=l 
77     TEST=Q** (K+l) / (P*DFLQAT (K+2) ) 
IF  (TEST.LE.EPS)  GO  TO  88 
K=K+1 
GO  TO  77 
88     M=K 

DETLAM=0.0D0 
DO  51  N=1,M 
RATIONAL. 0D0 
DO  52  JJ=1,N 
52     PJ^IC^PJ^ICW*DFLOAT(JJ)/(DFLQAT(JJ)+R) 
51     DETLAM=DETLAM+Q*  *N*RATION/DFLOAT  (N+l ) 
EFF2=-FCN  (TNOT) 
IF  (L.NE.l)  GO  TO  141 
UNOT=-DLOG (DABS (TNOT) ) 

WRTTE(6,98)  PSTAR,RSTAR,P/R/EFFl,EFF2,UNOT 
98     FORMAT(/, 11X, 57 C'&'hAlSX,  'VALUES  OBTAINED  AT  THE  END  OF  FIRST  IT 
1ERATION'  ,/,12X,  'PTHIL1 ,3X,  'RTHTL'  ,3X,  'PSTAR'  ,3X,  'RSTAR'  ,3X,  'ETHTL' 
2,3X, 'ESTAR' ,3X,  ,UNOT',/,12X,F5.3,2X,F7.3,2X,F5.3,2X,F7.3,2X,F5.3,3 
3X,F5.3,2X/F7.4///11X, 57 ('&'),/) 
141    IF  (DABS (P-PSTAR) .GT. 0.001)  GO  TO  29 
IF  (DABS (R-RSTAR)  .GT. 0.005)  GO  TO  29 
IF (DABS (EFF1-EFF2) .GT. 0.01)  GO  TO  29 
UNGT=-DLOG (DABS (TNOT) ) 

WRITE  (6,94)  XBAR,SSQR,P/R,EFF2,UNOT,L 
94     FORMAT  (15X,  'USING  THE  ABOVE  DATA  FOLLOWING  INFORMATION  IS  OBTAINED 
l///,15X,54('-'),/,10X,  'SPLE.MEAN',2X/'SPLE.VAR,,3X,  'EST. OF  P',3X,  ' 
2EST.OF  R'^X/EFF'^X,  'OP.UNOT' ,3X, 'NO.GF  ITR' ,/,10X,7  (9  ('-' )  ,2X)  , 
3/,10X,F9.3,2X,F9.3/2X,F9.3,2X,F9.3,2X,F9.3,2X,F9.4,5X,I3,//) 
GO  TO  114 
29     L=L+1 

GO  TO  95 

114    CONTINUE 

113    CONTINUE 

112    CONTINUE 

STOP 

END 
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SUBROUTINE  GNRNB(N,OBS,P,R,ISEED2,ISEEDl) 

KEAL*8  P,R,OBS 

DIMENSION  OBS (200) ,GAM(200) ,K(1) 

P1=SNGL(P) 

R1=SNGL(R) 

CALL  GAMA(R1,ISEED1,GAM,N) 

DO  1  1=1, N 

RIAM=GAM(I) * (1.-P1)/P1 

CALL  GGPOSH(RIAM,ISEED2,l,K/IER) 

OBS (I) =DBLE (FLOAT (K (1) ) ) 

RETURN 

END 


SUBROUTINE  COMPR(T,TX,XB,RST) 
IMPLICIT  REAL*8  (A-H,0-Z) 
COMMON  P,R,DETIAM 
R1=RST 
30     B=R1+XB* (1 . 0D0-T) 
TY=DABS (Rl/B) **Rl 
F=TY-TX 

FP1=TY 

FP2= (XB* (1 . 0D0-T) /B) +DLOG (DABS (Rl/B) ) 

FP=FP1*FP2 

R2=R1-F/FP 
IF(DABS(R2-R1) .LE.0.0001D0)  GO  TO  66 
R1=R2 
GO  TO  30 
66     R=R2 

RETURN 
END 
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